library(terra)
ndvi2018 = rast("data/sentinel2_ndvi2018s.tif")
ndvi2023 = rast("data/sentinel2_ndvi2023s.tif")
panel(c(ndvi2018, ndvi2023))Time for the above code chunk to run: 10.37 seconds
This document shows a set of examples related to the article “Comparing spatial patterns in raster data using R”. The examples are focused on comparing two layers of spatial continuous raster data, and are divided into three groups: raster outcome, single value outcome, and multiple value outcome.
The examples in this document use two raster layers representing NDVI values derived from Sentinel-2 images (“Sentinel-2 MSI Level-2A BOA Reflectance” 2022) for early summer of 2018 and 2023 in the area of Tartu, Estonia. The terra package (Hijmans 2024) is used for raster data manipulation.
The calculation time for each example is estimated and displayed at the end of each code chunk.
library(terra)
ndvi2018 = rast("data/sentinel2_ndvi2018s.tif")
ndvi2023 = rast("data/sentinel2_ndvi2023s.tif")
panel(c(ndvi2018, ndvi2023))Time for the above code chunk to run: 10.37 seconds
ndvi_diff = ndvi2023 - ndvi2018
plot(ndvi_diff)Time for the above code chunk to run: 1.23 seconds
ndvi_cor = focalPairs(c(ndvi2023, ndvi2018), w = 5,
fun = "pearson", na.rm = TRUE)
plot(ndvi_cor)Time for the above code chunk to run: 5.81 seconds
This example uses the geodiv package (Smith et al. 2023).
# 'sa': average surface roughness
library(geodiv)
window = matrix(1, nrow = 5, ncol = 5)
ndvi2018_sa_mw = focal_metrics(ndvi2018, window = window,
metric = "sa", progress = FALSE)
ndvi2023_sa_mw = focal_metrics(ndvi2023, window = window,
metric = "sa", progress = FALSE)
ndvi_diff_sa_mw = ndvi2023_sa_mw$sa - ndvi2018_sa_mw$sa
plot(ndvi_diff_sa_mw)Time for the above code chunk to run: 43.49 seconds
This example uses the GLCMTextures package (Ilich 2020).
library(GLCMTextures)
ndvi2018_q = quantize_raster(ndvi2018, n_levels = 16, method = "equal prob")
ndvi2023_q = quantize_raster(ndvi2023, n_levels = 16, method = "equal prob")
ndvi2018_textures = glcm_textures(ndvi2018_q, w = c(5, 5), na.rm = TRUE,
metrics = "glcm_homogeneity",
n_levels = 16, quantization = "none")
ndvi2023_textures = glcm_textures(ndvi2023_q, w = c(5, 5), na.rm = TRUE,
metrics = "glcm_homogeneity",
n_levels = 16, quantization = "none")
ndvi2023_textures_diff = ndvi2023_textures - ndvi2018_textures
plot(ndvi2023_textures_diff)Time for the above code chunk to run: 211.59 seconds
ndvi_diff = ndvi2023 - ndvi2018
ndvi_diff_autocor = autocor(ndvi_diff, method = "moran", global = FALSE)
plot(ndvi_diff_autocor)Time for the above code chunk to run: 1.02 seconds
This example uses the SSIMmap package (Ha and Long 2023).
library(SSIMmap)
ndvi_ssim = ssim_raster(ndvi2018, ndvi2023, global = FALSE, w = 5)
plot(ndvi_ssim)Time for the above code chunk to run: 5.5 seconds
This example uses the rasterdiv package (Rocchini et al. 2021).
library(rasterdiv)
ndvi2018_int = ndvi2018 * 100
ndvi2023_int = ndvi2023 * 100
ndvi2018_rao = paRao(ndvi2018_int, window = 5, progBar = FALSE)
ndvi2023_rao = paRao(ndvi2023_int, window = 5, progBar = FALSE)
ndvi_rao_diff = ndvi2023_rao[[1]][[1]] - ndvi2018_rao[[1]][[1]]
plot(ndvi_rao_diff)Time for the above code chunk to run: 556.63 seconds
This example uses the yardstick package (Kuhn, Vaughan, and Hvitfeldt 2024).
library(yardstick)
ndvi_rmse = rmse_vec(values(ndvi2023)[,1], values(ndvi2018)[,1])
ndvi_rmse[1] 0.2191853
Time for the above code chunk to run: 0.07 seconds
This example uses the diffeR package (Pontius Jr. and Santacruz 2023).
library(diffeR)
ndvi_mad = MAD(ndvi2023, ndvi2018)[1] "The mean of grid1 is less than the mean of grid2"
ndvi_mad Resolution Quantity Strata Element Total
1 10 0.1783152 0 0.005990721 0.184306
Time for the above code chunk to run: 0.36 seconds
This example uses the philentropy (HG 2018) package.
library(philentropy)
softmax = function(x) {
exp_x = exp(x - max(x))
return(exp_x / sum(exp_x))
}
ndvi2018_vals = na.omit(values(ndvi2018)[,1])
ndvi2023_vals = na.omit(values(ndvi2023)[,1])
ndvi2018_vals_prob = softmax(ndvi2018_vals)
ndvi2023_vals_prob = softmax(ndvi2023_vals)
ndvi2023_vals_prob_interp = approx(seq_along(ndvi2023_vals_prob),
ndvi2023_vals_prob,
xout = seq_along(ndvi2018_vals_prob))$y
ndvi_mat = rbind(ndvi2023_vals_prob_interp, ndvi2018_vals_prob)
philentropy::distance(ndvi_mat, method = "kullback-leibler")kullback-leibler
0.03759518
Time for the above code chunk to run: 0.11 seconds
This example uses the geodiv package (Smith et al. 2023).
library(geodiv)
ndvi2018_sa = sa(ndvi2018)
ndvi2023_sa = sa(ndvi2023)
abs(ndvi2023_sa - ndvi2018_sa)[1] 0.06886273
Time for the above code chunk to run: 0.05 seconds
This example uses the GLCMTextures package (Ilich 2020).
library(GLCMTextures)
ndvi2018_q = quantize_raster(ndvi2018, n_levels = 16, method = "equal prob")
ndvi2023_q = quantize_raster(ndvi2023, n_levels = 16, method = "equal prob")
ndvi2018_q_mat = as.matrix(ndvi2018_q, wide = TRUE)
ndvi2023_q_mat = as.matrix(ndvi2023_q, wide = TRUE)
ndvi2018_horizontal_glcm = make_glcm(ndvi2018_q_mat, n_levels = 16,
shift = c(1,0), na.rm = TRUE)
ndvi2023_horizontal_glcm = make_glcm(ndvi2023_q_mat, n_levels = 16,
shift = c(1,0), na.rm = TRUE)
ndvi2018_hom = glcm_metrics(ndvi2018_horizontal_glcm, "glcm_homogeneity")
ndvi2023_hom = glcm_metrics(ndvi2023_horizontal_glcm, "glcm_homogeneity")
abs(ndvi2023_hom - ndvi2018_hom)glcm_homogeneity
0.01880257
Time for the above code chunk to run: 0.5 seconds
This example uses the SSIMmap package (Ha and Long 2023).
library(SSIMmap)
ssim_raster(ndvi2018, ndvi2023, global = TRUE)SSIM: 0.63845 SIM: 0.90432 SIV: 0.90778 SIP: 0.75671
Time for the above code chunk to run: 0.85 seconds
This example uses the SpatialPack package (Vallejos, Osorio, and Bevilacqua 2020).
library(SpatialPack)
ndvi2018_mat = as.matrix(ndvi2018, wide = TRUE)
ndvi2023_mat = as.matrix(ndvi2023, wide = TRUE)
ndvi_CQ = CQ(ndvi2018_mat, ndvi2023_mat)
ndvi_CQ$CQTime for the above code chunk to run: 0 seconds
This example uses the bespatial package (Nowosad 2024).
library(bespatial)
ndvi2018_bes = bes_g_gao(ndvi2018, method = "hierarchy", relative = TRUE)
ndvi2023_bes = bes_g_gao(ndvi2023, method = "hierarchy", relative = TRUE)
abs(ndvi2023_bes$value - ndvi2018_bes$value)[1] 15402.11
Time for the above code chunk to run: 0.18 seconds
This example uses the keras3 (Kalinowski, Allaire, and Chollet 2024) and philentropy (HG 2018) packages.
library(keras3)
library(philentropy)
# keras3::install_keras(backend = "tensorflow")
normalize_raster = function(r) {
min_val = terra::global(r, "min", na.rm = TRUE)[[1]]
max_val = terra::global(r, "max", na.rm = TRUE)[[1]]
r = terra::app(r, fun = function(x) (x - min_val) / (max_val - min_val))
return(r)
}
ndvi_2018n = normalize_raster(ndvi2018)
ndvi_2023n = normalize_raster(ndvi2023)
ndvi2018_mat = as.matrix(ndvi_2018n, wide = TRUE)
ndvi2023_mat = as.matrix(ndvi_2023n, wide = TRUE)
ndvi2018_mat = array(rep(ndvi2018_mat, 3),
dim = c(nrow(ndvi2018_mat), ncol(ndvi2018_mat), 3))
ndvi2023_mat = array(rep(ndvi2023_mat, 3),
dim = c(nrow(ndvi2023_mat), ncol(ndvi2023_mat), 3))
model = keras3::application_vgg16(weights = "imagenet", include_top = FALSE,
input_shape = c(nrow(ndvi2018_mat), ncol(ndvi2018_mat), 3))
ndvi2018_mat = keras3::array_reshape(ndvi2018_mat, c(1, dim(ndvi2018_mat)))
ndvi2023_mat = keras3::array_reshape(ndvi2023_mat, c(1, dim(ndvi2023_mat)))
features2018 = predict(model, ndvi2018_mat)1/1 - 3s - 3s/step
features2023 = predict(model, ndvi2023_mat)1/1 - 2s - 2s/step
feature_map_2018_1 = as.vector(features2018[1,,,1])
feature_map_2023_1 = as.vector(features2023[1,,,1])
distance(rbind(feature_map_2018_1, feature_map_2023_1))euclidean
3.446775
Time for the above code chunk to run: 15.43 seconds
ndvi_diff = ndvi2023 - ndvi2018
hist(ndvi_diff)Time for the above code chunk to run: 0.15 seconds
This example uses the waywiser package (Mahoney 2023).
library(waywiser)
cell_sizes = c(200, 400, 600)
ndvi_multi_scale = ww_multi_scale(truth = ndvi2018, estimate = ndvi2023,
metrics = list(yardstick::rmse),
cellsize = cell_sizes,
progress = FALSE)
ndvi_multi_scale# A tibble: 3 × 6
.metric .estimator .estimate .grid_args .grid .notes
<chr> <chr> <dbl> <list> <list> <list>
1 rmse standard 0.193 <tibble [1 × 1]> <sf [625 × 5]> <tibble [0 × 2]>
2 rmse standard 0.191 <tibble [1 × 1]> <sf [169 × 5]> <tibble [0 × 2]>
3 rmse standard 0.188 <tibble [1 × 1]> <sf [81 × 5]> <tibble [0 × 2]>
Time for the above code chunk to run: 3.47 seconds